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t/3 Abstract 
O 

'55 Over the last few decades, climate scientists have devoted much effort to 

the development of large numerical models of the atmosphere and the ocean. 
While there is no question that such models provide important and useful in- 
formation on complicated aspects of atmosphere and ocean dynamics, skillful 
prediction also requires a phenomenological approach, particularly for very slow 
^ processes, such as glacial-interglacial cycles. Phenomenological models are often 

represented as low-order dynamical systems. These are tractable, and a rich 
source of insights about climate dynamics, but they also ignore large bodies 
of information on the climate system, and their parameters are generally not 
operationally defined. Consequently, if they are to be used to predict actual 
climate system behaviour, then we must take very careful account of the uncer- 
0^ tainty introduced by their limitations. In this paper we consider the problem of 

the timing of the next glacial inception, about which there is on-going debate. 
^ Our model is the three-dimensional stochastic system of Saltzman and Maasch 

^ (1991), and our inference takes place within a Bayesian framework that allows 

5-H both for the limitations of the model as a description of the propagation of 

the climate state vector, and for parametric uncertainty. Our inference takes 
the form of a data assimilation with unknown static parameters, which we per- 
form with a variant on a Sequential Monte Carlo technique ('particle filter'). 
Provisional results indicate peak glacial conditions in 60,000 years. 
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1 Introduction 



L 'analyse mathematique pent deduire des phenomenes generaux et simples 
Vexpression des lois de la Nature; mais Vapplication speciale de ces lois a des effets 
tres-composes exige une longue suite d' observations exactes. 

Joseph Fourier (1768 - 1830) 

This quote by Joseph Fourier appeared first in the 'Discours preUminaire' of the 
Analytical theory of heat [1] . With this sentence, Joseph Fourier expresses the need 
for an inductive approach to complex physical phenomena at the macroscopic scale. 
He repeated it at least once, to conclude his Memoire sur les temperatures du globe 
terrestre et des espaces planetaires [2] in which Fourier formulates what is known 
today as the greenhouse effect. Fourier confesses that "the question of Earth's tem- 
perature is one of the most important and difficult of all the Natural Philosophy" [2] 
and solving it was one central motivation for the theory of heat. Clearly, Fourier 
had fully perceived the complex character of the climate system. How, two centuries 
later, do we cope with climate's complexity? Which mathematical analysis is the 
most appropriate to get the best out of observations? In this paper we suggest that 
the most complicated model is not necessarily the most useful. Our understanding 
of the climate system should rather be judged by our ability to build physically 
consistent models that can be calibrated and validated on observations. 

We will illustrate our case with reference to a debate currently taking place in the 
circle of Quaternary climate scientists. The climate history of the past few million 
years is characterised by repeated transitions between 'cold' (glacial) and 'warm' 
(interglacial) climates. The first modern men were hunting mammoth during the 
last glacial era. This era culminated around 20,000 years ago [3j and then declined 
rapidly. By 9,000 years ago climate was close to the modern one. The current 
interglacial, called the Holocene, should now be coming to an end, when compared 
to previous interglacials, yet clearly it is not. The debate is about when to expect the 
next glacial inception, setting aside human activities, which may well have perturbed 
natural cycles. 

On one side. Professor Bill Ruddiman carefully inspected and compared palaeo- 
environmental information about the different interglacial periods. This comparison 
let him to conclude that glacial inception is largely overdue ^ |5j . According to him, 
the Holocene was not supposed to be this long, but the natural glacial inception 

^ Mathematical analysis allows you to deduce Nature's laws from general and simple phenomena; 
but applying these laws to highly composite effects requires a long series of exact observations. 
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process was stopped by an anthropogenic perturbation that began as early as 8,000 
years ago (rice plantations and land management by antique civilisations). On the 
other side, Professor Andre Berger and colleagues developed a mathematical model 
of the climate system, rated today as a 'model of intermediate complexity' |6l [7] 
to solve the dynamics of the atmosphere and ice sheets on a spatial grid of 19 x 5 
elements, with a reasonably extensive treatment of the shortwave and longwave 
radiative transfers in the atmosphere. Simulations with this model led Berger and 
Loutre to conclude that glacial inception is not due for another 50,000 years, as 
long as the CO2 atmospheric concentration stays above 220 ppmv [8]. Who is 
right? Perhaps both — Crucifix and Berger argued that the two statements are not 
strictly incompatible [9]. But perhaps neither. Both Ruddiman and Berger judge 
that it is possible to predict climate thousands of years ahead, but is this a realistic 
expectation? Michael Ghil wondered "what can we predict beyond one week, for how 
long and by what methods?" in a paper entitled Hilhert problems of the geosciences 
in the 21st century [TO]. This is the motivation for the present article. 

2 Steps towards a dynamical model of palaeoclimates 

2.1 An inductive approach to complex system modelling 

A system as complex as climate is organised at different levels : clouds, cloud sys- 
tems, synoptic waves, planetary waves, pluri-annual oscillations such as El-Niho, 
glacial-interglacial cycles, and so on. These patterns constitute information that is 
susceptible of being modelled and predicted. It is not our purpose to explain here 
how, in general, patterns emerge in complex systems ([TI] is an up-to-date reference 
on the subject) but it is useful to have a few notions in mind. Schematically, spatio- 
temporal structures are created by instabilities (necessarily fed by some source of 
energy) , and destroyed by relaxation processes (return to equilibrium) . The balance 
results in large-scale persistent patterns. In the atmosphere, local hydrodynamical 
instabilities can result in planetary waves, such as the ones responsible for dominant 
north-westerly winds in Canada and south-westerly winds in Europe. 

In most natural cases, the mechanisms of instability growth are so numerous and 
intricate that they would completely defeat reductionist modelling approaches for 
all but the shortest time intervals. The is certainly true for climate. For example, 
Saltzman repeatedly insisted |12| [T3] on the fact that neither current observations 
nor modelling of the present state of the atmosphere can possibly inform us of the 
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ice-sheet mass balances with sufficient accuracy to predict their evolution at the 
timescale of several thousands of years. But, at the same time, the palaeoclimate 
archive suggests the behaviour of climate to be fairly regular, over hundreds of 
thousands of years. Consequently, if our interest is in this large-scale behaviour 
over long time intervals, then we may expect to do better by modelling it directly. 
An extreme case of this is simply to fit time-series models. However, these require 
detailed and accurate records, and have only a limited ability to account for the 
influence of time-varying factors, like, in the case of the climate, solar insolation 
arising from orbital variations. Between the two extremes of a reductionist model 
and a time series model there is a huge space for physically-based models that 
may be empirically tuned. These are both physical models, in that they embody 
basic principles in the way that they constrain the interaction between variables, and 
statistical models, since the coefficients in these models may be treated as uncertain, 
to be learnt about from observations. In physics, models that treat observable 
quantities directly are sometimes termed 'phenomenological' . Our intention is to 
go further, and embed the phenomenological model in a statistical framework. This 
allows us to replace the somewhat haphazard tuning process that these models often 
undergo with a probability-based updating process, as we will describe in section 3. 

2.2 Empirical evidence about the Quaternary 

Building a phenomenological model of glacial-interglacial cycles requires a detailed 
knowledge of the Quaternary. This section provides a brief overview of the vast 
amount of information that scientists have accumulated on that period before we 
tackle Ruddiman's hypothesis. 

2.2.1 The natural archives 

By the 1920s, geomorphologists were able to interpret correctly the glacial moraines 
and alluvial terraces as the left-overs of previous glacial inceptions. Penk and 
Briickner ([H], cited by [15]) recognised four previous glacial epochs, named the 
Giinz, Mindel, Riss and Wiirm, taking us back about 400,000 years. The wealth of 
data on the Quaternary environments that has since been collected and analysed 
by field scientists can be appreciated from the impressive four-volume Encyclopedia 
of Quaternary Sciences recently edited by Elias flG]. Analysing and interpreting 
palaeoenvironmental data involves a huge variety of scientific disciplines, including 
geochemistry, vulcanology, palaeobiology, nuclear physics, stratigraphy, sedimentol- 
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ogy, glacial geology and ice-stream modelling. Here is a brief summary. 

Stable isotopes constitute one important class of natural archives. It has 
been known since the works of Urey [T7], Buchanan [18] and Dansgaard [19] that 
physical and chemical transformations involved in the cycles of water and carbon 
fractionate the isotopic composition of these elements. For example: ice-sheet water 
is depleted in oxygen-18 and deuterium compared to sea water; clouds formed at low 
temperatures are more depleted in oxygen-18 and deuterium than clouds formed at 
higher temperatures; organic matter is depleted in ^'^C, such that inorganic carbon 
present in biologically active seas and soils is enriched in ^'^C (^^N is another useful 
stable palaeo-environmental indicator sensitive to the biological activity of soils). 
The isotopic compositions of water and biogenic carbon are extracted from deep- 
sea sediments, ice and air trapped in ice bubbles, palaeosols, lake-sediments and 
speleothemes. One of the first continuous deep-sea record of glacial- interglacial 
cycles was published by Cesare Emiliani |20j . 

Radioactive tracers are used to estimate the age of the record and rate of 
ocean water renewal. At the timescale of the Quaternary, useful mother-daughter 
pairs are ^sorpj^ j 238, 234 ||{ja,ting carbonates), and ^''K / ^"^Ar in potassium-bearing 
minerals. The ratio ^SOrpj^ j 23ipg^ jg ^ useful indicator of ocean circulation rates. 

The chemical composition of fossils is also indicative of past environmental 
conditions. In the ocean, the amounts of cadmium, lithium, barium and zinc trapped 
in the calcite shells of foraminifera indicate the amount of nutrients at the time of 
calcite formation, while the amounts of magnesium and strontium are empirically 
related to water-temperature. 

Glaciologists have also developed ambitious programmes to analyse the com- 
position of air (oxygen, nitrogen, plus trace gases such as methane, carbon-dioxide 
and nitrogen oxide, argon and xenon) trapped in ice accumulating on ice sheets, of 
which the European Project for lee Core in Antarctica (EPICA) is a particularly no- 
table achievement [21j. It was demonstrated that the central plateaus of Antarctica 
offer a sufficiently stable environment to preserve air's chemical composition over 
several hundreds of thousands of years. The longest presently available trace gas 
record — including carbon dioxide — goes back to 800 kyr [22] 

Other sources of information. Plant and animal fossils (including pollens) 
trapped in lakes, peat-bogs, palaeosols and marine sediments indicate the then- 
prevailing palaeoenvironmental conditions. Their presence or absence may be inter- 
preted quantitatively to produce palaeoclimatic maps [23|. Preservation indicators of 
ocean calcite fossils are used to reconstruct the history of ocean alkalinity. Palaeosols 
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and wind-blown sediments (loess) indicate past aridity. The loess grain-size distribu- 
tion is also sensitive to atmospheric circulation patterns. Geomorphological records 
are an important source of information about the configuration of past ice sheets, 
which is complemented by datable evidence (typically coral fossils) on sea-level. 

2.2.2 The structure of Quaternary climate changes 

It is not straightforward to appreciate which aspects of a climate record are relevant 
to understand climate dynamics at the global scale. For example, minor shifts in 
oceanic currents may affect the local isotopic composition of water, but with no 
observable effects on glacial-interglacial cycle dynamics. One strategy is to collect 
samples from many areas of the world and average them out according to a process 
called 'stacking'. One of the first stacks, still used today, was published by John 
Imbrie and colleagues [23] in the framework of the Mapping spectral variability in 
global climate project; it is usually referred to as the Specmap stack. Here we con- 
centrate on the more recent compilation provided by Lisiecki and Raymo [26] , called 
LR04. The stack was obtained by superimposing 57 records of the oxygen-18 com- 
position of benthic foraminifera shells. Benthic foraminifera live in the deep ocean 
and therefore record the isotopic composition of deep water (an indicator of past ice 
volume). However, there is an additional fractionation associated with the calcifi- 
cation process, which is related to water temperature. The isotopic composition of 
calcite oxygen is reported by a value, named 5^^0c, giving the relative enrichment 
of oxygen-18 versus oxygen-16 in calcite (hence the 'c' subscript) compared to an 
international standard. High 6^^0 indicates low continental ice volume and/or high 
water-temperature. 

The LR04 stack (Figure [T]) shows the gradual transition from the Pliocene — 
warm and fairly stable — to the spectacular oscillations of the late Pleistocene. The 
globally averaged temperature at the early Pliocene was about 5°C higher than today 
([25] and references therein); that at the last glacial maximum (20,000 years ago) 
was about 5°C lower. Key research questions are to characterise these oscillations, 
understand their origin and quantify their predictability. 

The Morlet Continuous wavelet transform |i28j provides a simple summary of the 
backbone of these oscillations. The LR04 record is dominated for most of the time 
by a 40,000 yr signal until roughly 900,000 years ago, after which the 40,000 yr signal 
is still present but topped by longer cycles. At the very least, this picture demon- 
strates that LR04 contains structured information susceptible of being modelled and 
possibly predicted. 
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Figure 1: The LR04 benthic S^^O stack constructed by the graphic correlation of 57 globally 
distributed benthic S^^O records Note that the full stack goes back in time to —5.2 Myr 
(1 Myr = 1 million years, lkyr=1000 years). The signal is a combination of global ice volume 
(low S^^O corresponding to low ice volume) and water temperature (low S^^O corresponding 
to high temperature). Top: The time series, as available from www. loraine-lisiecki . com. 
The vertical axis is reversed as standard practice to get 'cold' climates down. Bottom: 
Modulus of its continuous Morlet wavelet transform (algorithm of [55]) using luq = 5.4 but 
with a normalisation c(s) = s/y/At. Large amplitudes are in black. 
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How complex a model is needed to model such a series? There is no clear-cut 
answer. Time-series extracted from complex systems are sometimes characterised 
by their correlation dimension, which is an estimator for the fractal dimension of the 
corresponding attractor [30 . The first estimates for the Pleistocene were provided 
by Nicohs and NicoHs ^ (d = 3.4) and Maasch et al. [32] (4 < d < 6). For this 
article we calculated correlation dimension estimates for the LR04 stack (d ~ 2.4) 
and the HW04 stack [33j {d ~ 2.5); estimates were made using the 'fractal' R 
package [3l] . HW04 is similar to LR04 but it is based on different records and dating 
assumptions. Note, however, that several authors, including Grassberger |35 l l36 l 157] . 
have discouraged the use of correlation dimension estimates for the 'noisy and short' 
time series typical of the Quaternary because they are overly sensitive to sampling 
and record length. 

In response to this problem Ghil and colleagues [ST] [38] advocated singular- 
spectrum analysis (SSA), in which a time series is linearly decomposed into a number 
of prominent modes (which need not be harmonic), plus a number of small-amplitude 
modes. Assuming that the two groups are separated by an amplitude gap, the first 
group provides the low-order backbone of the signal dynamics while the second group 
is interpreted as noise. SSA was applied with some success to various sediment and 
ice-core records of the last few glacial-interglacial cycles [38] and has in general 
confirmed that the backbone of climate oscillations may be captured as a linear 
combination of a small number of amplitude and/or frequency-modulated oscilla- 
tions, plus a long-term trend. SSA of the last 900 kyr of LR04 (Figure |2]) confirms 
this statement, but ordering eigenvalues according to their characteristic period in- 
dicates that the 'noisy background' of LR04 is essentially organised as first-order 
autoregressive process. 

2.2.3 The Achilles heel of reconstructions 

Now we must mention a particularly difficult and intricate issue: dating uncertainty 
in palaeoclimate records. No palaeoclimate record is dated with absolute confidence. 
Marine sediments are coarsely dated by identification of a number of reversals of 
Earth's magnetic field, which have been previously dated in rocks by radiometric 
means ([27^ and references therein). Magnetic reversals are rare (four of them over 
the last 3 million years) and their age is known with a precision no better than 5,000 
years. Local sedimentation rates may vary considerably between these time markers, 
such that any individual event observed in any core taken in isolation is hard to date. 
Irregularities in the sedimentation rate blur information in the spectral analysis. 
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Figure 2: Singular Spectrum Analy- 
sis (SSA) of the LR04 benthic stack. 
The figure shows the eigenvalues of 
the lagged-covariance matrix of rank 
M = 100 as given by [35], eq. (6). 
The record was cubic-spline interpo- 
lated {At = Ikyr) and only the most 
recent 900 kyr were kept. Eigen- 
values are here plotted according to 
their characteristic period. They 
then appear to be broadly consistent 
with a Ist-order autoregressive pro- 
cess (the ARl-nuU hypothesis is rep- 
resented by the grey shade 001111]), 
but significant features do appear at 
periods of 20 kyr and longer. Note 
that the spectrum is conditioned by 
astronomical tuning assumptions, as 
may be observed by comparison with 
non-tuned records. 



One strategy to tackle this issue is to assume synchrony between oscillation pat- 
terns identified in different cores. Statistical tests may then be developed on the 
basis that dating errors of the different cores are independent. For example, Huy- 
bers (2007) [42^ considered the null-hypothesis that glacial-to-interglacial transitions 
(they are called terminations in the jargon of palaeoclimatologists) are independent 
of the phase of Earth's obliquity. While this null-hypothesis could not be rejected on 
the basis of a single record, the combination of 14 cores allowed him to reject it with 
% confidence, proving once more the effect of the astronomical forcing on climate 



(see section 2.3). Note statistical tests to assess the significance of a correlation 



between two ill-dated palaeoclimate records are only now being developed 

Another strategy is known as orbital tuning. The method consists in squeezing 
or stretching the time-axis of the record to match the evolution of one or a combi- 
nation of orbital elements, possibly pre-filtered by a climate model |24l HI] . This 
strategy has undeniably engendered important and useful results (e.g. [IS]), but at 
the cost of imposing link between orbital forcing and the record. Experienced inves- 
tigators recognise that orbital tuning has somehow contaminated most of the dated 
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palaeoclimate records available in public databases. Quantification of leads and lags 
between CO2 and ice volume in particular is hostage to hidden dating assumptions 
and circular reasoning. Here is one typical illustration. Saltzman and Verbitsky 
showed on several occasions (e.g. [46j) a phase diagram showing the SPECMAP 
6^^0 stack versus the first full ice-core records of CO2 from Vostok j471 148j . It is 
reproduced here (Figure |3j left). The phase diagram clearly suggests that CO2 leads 
ice volume at the 100 kyr time scale. However, a detailed inspection of the original 
publications reveals that the SPECMAP record was astronomically tuned, and that 
the Vostok time-scale, called GT4, uses a conventional date of isotopic stage 5.4 of 
110 kyr BP ... by reference to SPECMAP! ^8\. The observed hysteresis may there- 
fore be an artifact. The situation today is that there is no clear consensus about the 
phase relationship between ice volume and CO2 at the glacial-interglacial time scale 
(compare [121 [501 [51] ) . According to Ruddiman's analysis [50], CO2 leads ice volume 
at the precession (20 kyr) period, but CO2 and ice volume are roughly synchronous 
at the obliquity (40 kyr) period. Current evidence about the latest termination is 
that the decrease in ice volume and the rise in CO2 began around 19,000 years ago 
[I9l[52]. 

2.3 Getting physical laws into the model 

So far we have learned that palaeoclimate oscillations are structured and that it is 
not unreasonable to attempt to model them with a phenomenological model forced 
by the astronomical variations of Earth's orbit. What is the nature of the physical 
principles to be embedded in such a phenomenological model, and how can they 
be formalised? The history of Quaternary modelling is particularly enlightening in 
this respect ( [15] provides an extensively documented review of Quaternary climates 
modelling up to the mid-1980s). 

Joseph Adhemar (1797-1862) [5l] first suggested that the cause of glaciations is 
the precession of the equinoxes, but subsequently there was disagreement about how 
precession affected climate. Joseph John Murphy maintained that cold summers 
(occurring when summer is at aphelion) favour glaciation [55], while James Croll 
(1821-1890) maintained that cold winters are critical [56]. Croll's book demonstrates 
a phenomenal encyclopaedic knowledge. His judgements were in places particularly 
far-sighted, but they are barely substantiated by the mathematical analysis Fourier 
was so insistent about. The nature of his arguments was essentially empirical, if not 
in places frankly rhetorical. 
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Figure 3: The concentration in CO2 measured in the Vostok ice core record [53] over the last 
200 kyr is plotted versus two proxies of continental ice volume. Left: the planctonic S^^O 
stack by Imbrie et al. (1984) (standard units). Right: the benthic S^^O stack by Lisiecki 
and Raymo (2004) (per mil). Numbers are dates, expressed in kyr BP (before present). The 
Imbrie stack suggests hysteresis with CO2 leading ice-volume, but LR04 does not show this 
so clearly. 

Milutin Milankovitch (1879-1958) is today acknowledged as the pioneer of math- 
ematical climatology. In a highly technical book that crowned a series of articles 
written between 1920 and 1941 |57| . Milankovitch extended Fourier's work to esti- 
mate the zonal distribution of Earth's temperature from incoming solar radiation. 
He also computed the effects of changes in precession, eccentricity and obliquity 
on incoming solar radiation at different latitudes to conclude, based on geological 
evidence, that summer insolation drives glacial-interglacial cycles, thus supporting 
Murphy. 

Mathematical analysis allowed Milankovitch to deduce the consequences of cer- 
tain fundamental principles, such as the laws of Beer, Kirchhoff and Stefan, on global 
quantities such as the Earth's temperature. But Milankovitch also used empirical 
macroscopic information, such as the present-day distribution of the snow-line al- 
titude versus latitude, to estimate the effects of temperature changes on the snow 
cover. In today's language, one may say that Milankovitch had accepted that some 
information could not be immediately inferred from microscopic principles because it 
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depended on the way the system as a whole responded to its numerous and intricate 
constraints (Earth's rotation, topography, air composition etc.). 

The marine-record study pubhshed by Hays, Imbrie and Shackleton |44] is often 
cited as the most indisputable proof of Milankovitch's theory. Hays et al. identified 
three peaks in the spectral estimate of climate variations that precisely correspond 
to the periods of obliquity (40 kyr) and precession (23 kyr and 19 kyr) calculated 
analytically by Andre Bergeij^ Milankovitch's theory of ice ages, though, had to be 
revisited in the face of evidence — already available in an article by Broecker and van 
Donck [62] — that the glacial cycle is 100,000 years long, ice build-up taking about 
80,000 years and termination about 20,000 years |111|62]. Neither the 100,000-year 
duration of ice ages, nor their saw-tooth shape were predicted by Milankovitch. The 
bit Milankovitch's theory is missing is the dynamical aspect of climate's response. 
Glaciologist Johannes Weertman [63] consequently addressed the evolution of ice 
sheet size and volume by means of an ordinary differential equation, thereby open- 
ing the door to the use of dynamical system theory for understanding Quaternary 
oscillations. 

In the meantime, general circulation models (GCMs) of the atmosphere and 
oceans running on supercomputers became widely available (cf. [6^ for a review), 
and used for palaeoclimate purposes |65 | I66 | [67]. Just as Milankovitch applied Beer 
and Kirchoff's laws to infer Earth's temperature distribution, GCMs allow us to 
deduce certain aspects of the global circulation from our knowledge of balance equa- 
tions (i.e. representing continuity and conservation) in each grid cell. However, these 
balance equations are highly uncertain at the resolution of current solvers — clouds, 
for example, are very hard to represent reliably — and quantifying this uncertainty 
at the global scale is a difficult problem which is only now being systematically 
addressed j68t [BH] . While GCMs are undeniably useful to constrain the immediate 
atmospheric response to changes in orbital parameters, they are far too uncertain 
to estimate glacial accumulation rates reliably enough to predict the evolution of ice 
sheets over tens of thousands of years [12] . 

In the following sections we will concentrate on a three-dimensional dynamical 
climate model proposed by Barry Saltzman. This choice was guided by the ease of 
implementation as well as the impressive amount of supporting documentation |13) . 
However, there were alternatives to this choice. The reader is referred to the article 

■^The supporting papers by Berger would only appear in the two following years [5811591 15D] . Hays 
et al. based their analysis on a numerical spectrum estimate of the orbital time-series provided by 
Vernekar [61| . 
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by Imbrie et al. [70] and pp. 264-265 of Saltzman's book [13] for a summary with 
numerous references organised around the dynamical concepts proposed to explain 
glacial-interglacial cycles (linear models, with or without self-sustained oscillations, 
stochastic resonance, model with large numbers of degrees of freedom). 

Among the alternatives to Saltzman's model, the series of models published by 
Ghil and colleagues [TTJ [721 E3] have some of the richest dynamics. They present 
self-sustained oscillations with a relatively short period (6,000 years). The effects 
of the orbital forcing are taken into account by means of a multiplicative coefficient 
in the ice mass balance equation. This causes non-linear resonance between the 
model dynamics and the orbital forcing. The resulting spectral response presents 
a rich background with multiple harmonics and band-limited chaos. More recently, 
Gildor and Tziperman [71] proposed a model where sea-ice cover plays a central 
role. In this model, termination occurs when extensive sea-ice cover reduces ice 
accumulation over ice sheets. Like Saltzman's, this model presents 100 kyr self- 
sustained oscillations that can be phase-locked to the orbital forcing. We believe 
that much is to be learned from a systematic analysis of these simple models, and 
they easily find their place in the climate modelling hierarchy proposed by Held [75] . 

Scientists with long field experience have also proposed models. These are usually 
qualified as 'conceptual', in the sense that they are formulated as a worded causal 
chain inferred from a detailed inspection of palaeoclimate data without the support 
of differential equations. Good examples are [50] (TO] [TH [77]. In the two latter 
references, Ruddiman proposes a direct effect of precession on CO2 concentration and 
tropical and southern-hemisphere sea-surface temperatures, while obliquity mainly 
affects the hydrological cycle and the mass-balance of northern ice sheets. 

3 The Saltzman model (SM91) 
3.1 Outline 

As a student of Edward Lorenz, Barry Saltzman (died 2002) contributed to the 
formulation and study of the famous Lorenz63 dynamical system [TH] traditionally 
quoted as the archetypal low-order chaotic systerrj^ Saltzman was therefore in an 
excellent position to appreciate the explanatory power of dynamical system theory. 

^The Acknowledgments of the Lorenz (1963) paper read: "The writer is indebted to Dr. Barry 
Saltzman for bringing to his attention the existence of nonperiodic solutions of the convection 
equations." 
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Between 1982 and 2002 he and his students pubHshed several dynamical systems 
deemed to capture and explain the dynamics of Quaternary oscillations |12| [T3l [79l 
[Hoi ED US] • In the present article we choose to analyse the 'palaeoclimate dynamical 
model' published by Saltzman and Maasch (1991) [81]. We will refer to this model 
as SM91. 

Saltzman judged that the essence of Quaternary dynamics should be captured 
by a three-degree-of-freedom dynamical system, possibly forced by the variations in 
insolation caused by changes in orbit [12]. The evolution of climate at these time 
scales is therefore represented by a trajectory in a 3-dimensional manifold, which 
Saltzman called the "central manifold". The three variables are ice volume (/), 
atmospheric CO2 concentration and deep-ocean temperature (9). It is important 
to realise that Saltzman did not ignore the existence of climate dynamics at shorter 
and longer time scales than those that characterise the central manifold, but he 
asserted the hypothesis that these modes of variability may be represented by distinct 
dynamical systems. In this approach, the fast relaxing modes of the complex climate 
system are in statistical equilibrium with its slow and unstable dynamical modes. 
This is sometimes called the 'adiabatic elimination' |83] . 

A statistical representation of the fast modes is plainly justified when a spectral 
gap separates slow and fast variances. Unfortunately — and contrary to a common 
perception |84| l85] — there is no evidence for such a gap in the palaeoclimate record 
\86\ IHTl |88] . However, the latter references document a change in the power-law 
exponent of the spectral background: signal energy decays faster with frequency 
below (lOOyr)"^ than above. This is evidence that the effective dissipation time- 
scale is larger for high frequencies, and that the dynamics of slow and fast climatic 
oscillations are at least partly decoupled. To be complete, note that [86] contrasted 
atmospheric variability during modern glacial era on the basis Greenland ice core 
data, and found that atmospheric variance was larger, and fatter-tailed during the 
glacial era. 



3.2 Formulation 

The ice mass-balance is the result of the contribution of four terms: a drift, a term 
inversely proportional to the deviation of the summer mean temperature at high 
latitudes (r) and a relaxation term 

^ = (^1 - - ^sl- (1) 
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By diabatic elimination, r is in thermal equilibrium with the slow variables [/, fi, 9] 
and its mean may therefore be estimated as a diagnostic function of the latter. 
Saltzman adopted a linear approximation of the form: 

T = To + ai{I - Iq) + a^{n - Ho) + ae{0 - Oq) + grR, (2) 

where R is the astronomical forcing. Saltzman used incoming insolation at 65° N 
at summer solstice, which can be calculated using the Berger algorithm |60] . The 
reference state [Iq, fiQ, 6q] will be defined later on. 

The CO2 equation includes the effects of ocean outgassing as temperature in- 
creases, a forcing term representing the net balance of CO2 incorporated into the 
atmosphere minus that eliminated by silicate weathering, and a non-linear dissipa- 
tive term: 

^ = f3o-l3ee + F^-K^fi with = - p^,^ + (3^^^. (3) 

The dissipative term (K^fi) is a so-called Landau form and its injection into the CO2 
equation is intended to cause instability in the system. In earlier papers (e.g. [89]), 
Saltzman [and Maasch] attempted to justify similar forms for the CO2 equation on a 
reductionist basis: each term of the equation was identified to specific, quantifiable 
mechanisms like the effect of sea-ice cover on the exchanges of CO2 between the 
ocean and the atmosphere, or the ocean circulation on nutrient pumping. It is 
noteworthy that Saltzman gradually dropped and added terms to this equation 
(compare |80| l8H [89] ) to arrive at the above formulation in which they posit a 
carbon cycle instability without enquiring too deeply into causal mechanisms. 

The deep-ocean temperature simply assumes a negative dependency on ice vol- 
ume with a dissipative relaxation term: 

d9 

= 71 - 72I - (4) 

Without loss of generality we are now free to assume that [Iq^iiq^Oq] is a fixed 
point of the evolution equations. Then evolution equations for [/', /i', 6'\ = [/, /z, 9\ — 
[Io,/io,^o] may be expressed as follows: 
dl' 

— = -ai [k^fi' + keO' + kRR'it)] - Kjl' (5a) 
^=6i/x'-62/x'2-63^'3-6e0' (5b) 

JO/ 

^ = -c^I' - K',e' (5c) 



15 



The different parameters ai,bi,ci,ki are functions of the Oi, fii,^i, Iq, fiQ and Oq. 
Some of these we will treat as known and fixed, but others we will treat as unknown, 
as available to be tuned to observations. For this we need a formal statistical frame- 
work that unifies the SM91 model and the climate system that the model purposes 
to represent. 

4 The Bayesian inference process 

Figure |4] reproduces the original solution [81] , using the parameters published at the 
time. As in the original publication, the solution is compared with Imbrie's (5^^0- 
stack |2l] interpreted as a proxy for ice volume, but we have added here the CO2 
record extracted from the Vostok and EPICA (Antarctica) ice cores \22\ |SH] . Clearly 
the agreement between the model and the observations is seductive, but perhaps we 
should not be surprised by this, given that, many poorer models might have been 
rejected, and that the parameters in this model have been tuned in some manner to 
optimise the fit. Several authors have also argued that it is difficult to chose among 
models with similar dynamical behaviours but built on different interpretations of 
the climate system's functioning |9T1 192] . So are such models useful ? 

Let us first observe that while this question is a pressing one for phenomenological 
models, it also concerns large general circulation models models (GCMs). Once 
assembled, GCMs are 'tuned' to capture major characteristics of climate such as 
the overturning cell or the global mean temperature (e.g. [93]). This tuning is an 
effective way of incorporating macroscopic information in the model, but then this 
information can no longer be treated as an emergent feature, i.e. an opportunity for 
validation, and a reason for confidence in prediction. In philosophy, this conundrum 
is sometimes termed the problem of prediction versus accommodation |94| [95] . 

We believe that the difficulties of validation and prediction can be mitigated 
by the use of a formal framework in which all choices are clearly specified, and 
which proceeds according to well-established rules. We will use a Bayesian statistical 
approach for this purpose. This approach has its roots in the works of Bayes, Laplace 
and Bernouilli, who were looking for ways to augment their knowledge of uncertain 
quantities such as initial conditions or parameters, by means of observations |96| . 
Rougier (2007) |97] explains how Bayesian inference methods may be applied to 
the problem of climate prediction. We reprise this explanation here, in a slightly 
more general form, and adapted to the problems posed by palaeoclimate time-series 
analysis. In particular, we explicitly consider the dynamical aspects of climate, the 
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Figure 4: Response of the palaeoclimate model of Saltzman and Maasch (1991) |81_. Shown 
are the insolation forcing, taken as the summer solstice incoming solar radiation at 65° N 
[BO] : the ice volume anomaly (full), overlain with the SPECMAP planctonic S^^Oc stack [53] 
(dashed), the CO2 atmospheric concentration, overlain with the Antarctic ice core data from 
Vostok and EPICA [221 [53], and finally deep-ocean temperature. SPECMAP foraminifera 
data are scaled linearly assuming that the transition between the last glacial maximum and 
today represents a melt-down of 44 • 10^^ kg of ice [3l|9D] (this scaling is less flattering as the 
one used by Saltzman, but it conforms better to the literature). Note that ice volume and 
ocean temperature are anomalies to the long-term (tectonic) average, and fi is calculated 
assuming ^0 = 253 ppmv. A similar figure was shown in the original article by Saltzman 
and Maasch [81]. 
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evolution of which is described by time-differential equations. 

Before embarking on the mathematical details, it is useful to recall two aspects in- 
herent in complex system modelling, as introduced in section 2.1 First, by focussing 
on certain modes of climate variability we ignore a large body of information, such 
as its synoptic variability, or, e.g.: the occurrence of a particular volcanic eruption 
at a particular time. This causes prediction errors that have to be incorporated 
stochastically, since they are uncertain (we will here neglect the epistemological dis- 
tinction between stochastic forcing and error). Thus validation consists in verifying 
that the complete set of model assumptions is compatible with observations, that is 
to say both the structure of physical model, and the way in which we incorporate its 
limitations. In general it is not possible to test either of these two parts in isolation. 
Second, our physical model is constructed using parameters that cannot in practice 
be deduced from our knowledge of microscopic interactions. In our case, these occur 
in conjectures about the mathematical expressions of carbon, ocean and ice-sheet 
feedbacks, which are to be calibrated, or tuned, by reference to observations. Thus 
we must treat these parameters as uncertain, in order that we can learn about them. 

Prior to the present study Hargreaves and Annan |[98] presented a Bayesian cali- 
bration of a slightly different version of the Saltzman and Maasch model. Our treat- 
ment addresses an acknowledged limitation of their approach, which is that they 
treat their model as perfect, except for uncertainty about the model-parameters. 
This is seen in their likelihood function, in which they compare the model trajectory 
at a given choice of parameters with the observed trajectory, and scale by observa- 
tional error (which is assumed to be independent across observations). We explicitly 
include a term in the model to account for the limitations of the model in describing 
the climate state at time t + dt given the climate state at time t. Unlike Hargreaves 
and Annan, we will not be able to give a simple closed-form expression for the like- 
lihood function of the model parameters, because once the discrepancy is included 
in this way, the joint distribution of the parameters and the evolving state vector 
has a very complicated structure, reflecting the underlying model. 

Now we outline the mathematical treatment. Denote by X{t) the state of the 
climate at time t, and hy X = {X{t)} the climate process over the period of interest. 
We do not observe X directly, but some proxy for it, which we denote as Z, and 
we denote the value of the observations as z. Typically the observations in Z are 
irregularly-spaced in time, and the relationship between X and Z is many-to-one; 
for example, Z might respond to a subset of the components of the state vector, or 
to a subset of spatial locations. 
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Our inference about X based on the observations Z = z can be expressed in 
probabilistic terms as 'k{X \ Z = z) where 7r(-) denotes a density function identified 
by its arguments, and '|' denotes 'conditional upon'. In other words, we would like 
to infer the density function of X conditional upon the fact that Z = z. Using the 
rules of the probability calculus, we can represent this inference in a way that is 
more amenable to calculations, 

7r(X I Z = z) oc 7r(Z = z I X)7r(X), (6) 

by applying Bayes's Theorem — the missing constant is the reciprocal of 7r(Z = z). 
The density 7r(Z = z \ X) describes the forward relationship between the climate 
state vector and the measured proxy: we will come back to this in section [5} 

The density t^{X) in ^ describes our uncertainty about X, and we introduce 
a physical model to help assess this density. In our case, we have a deterministic 
model for the evolution of X{t)^ which is specified conditional upon knowledge of 
the parameters -0. This model is written 

^=/(X(t),t,V^). (7) 

However, we are aware that this model has many limitations: some we know about, 
and could include but choose not to, others we do not even know. To capture the 
gross effect of these limitations, we introduce a stochastic component, writing our 
model in the Ito form 

dX{t) = f{X{t),t, V') dt + S^/^ j^(^) 

where S is a variance matrix that we must specify, and W{t) is a vector of inde- 
pendent Brownian motions. X \ ^ is now a stochastic process, and consequently 
X is a more complicated stochastic process, allowing ip to be uncertain. However, 
both sources of uncertainty are necessary, because even if we knew ip, we would not 
expect our model to do a perfect job of replicating the true system values X; and 
in fact we do not know ^p. Therefore to account for the limitations of our model we 
will need to specify a density function n^ip) and a variance matrix S. While neither 
of these tasks is easy, it is better to do the best we can, than to accept obviously 
wrong choices, such as asserting that ip is known, i.e. Tr^ip) = — tpo) for some ipQ 
that we specify, or that the model is perfect modulo our uncertainty about ip, i.e. 
setting S = 0. 
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Once we introduce our physical model we introduce the new set of uncertain 
quantities, ip. Hence ^ becomes 

tt{X, tp \ Z = z) (XTr{Z = z \ X) tt{X, ip) 

= Tr{Z = z \ X)tt{X \ i^)TT{%b), (9) 

where in the first line we assert that X is sufficient for Z, i.e. ^ is not relevant to Z 
when we have X, and the rearrangement in the second line is a standard probabilistic 
relationship. The second term in the second line is the density function of the 
stochastic model, (|8]), while the third term is our prior assessment of uncertainty 
about -0. In practice, we will not be able to evaluate the density function of ([s]), but 
we will be able to sample from it, which is sufficient for a numerical approximation 
of (§. 

Eq. Q is completely general: it requires no special conditions on the physical 
model, or on our probability assessments. However, our dynamical stochastic model 
does have special properties that we can exploit, notably that X \ is a Markov 
process, so that 

7r{X{t + l) I X{t),X{t-l),...,i,) =7T{X{t + l) I X{t),i;); (10) 

this can be summarised as stating that the past is irrelevant to the future, given the 
present. Miller |99] describes it thus: "Markov processes: the ordinary differential 
equations of the stochastic process realm" (p. 18). This dramatically simplifies the 
task of sampling from vr (X, ip \ Z = z), because it allows us to proceed sequentially 
through time. The Markov structure can be visualised using a Directed Acyclic 
Graph (DAG), as shown in Figure |4j In this DAG, the middle row is expressed 
by plus, if required, the initial density function 7r(X(0)): we will not require 
7r(X(0)) because we will spin- up the dynamical model from its long-term (hereafter 
termed 'tectonic') mean. The bottom row is expressed by 7r{ip). The connections 
between the middle and top rows express the forward relationship from the climate 
state vector at time t to the proxy measurement at time t, and we have treated the 
climate state at time t as sufficient for the proxy measurement at time t. 

Finally, we turn to the difficult issue of choosing ^T{^p) and S. As we have already 
stressed, the challenge here is not to specify the 'right' choices, but, at least initially, 
to avoid poor choices. We regard both ir^ip) = 5{tp — tpo) and S = as poor choices. 
It is very hard to make 'not poor' choices directly, because neither not S have 
clear operational definitions. We reject out-of-hand the notion that any choices can 
represent 'pure ignorance'. A better approach is to make choices on the basis of their 
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X(tp - 2) ^ - 1) ^ Xitp) ^ X{tp + 1) 




Figure 5: Directed Acyclic Graph of the inference, where tp is the last time-point for which 
we have proxy measurements. The process X(t) \ ip is Markov, allowing us to process the 
measurements sequentially in time. 

implications for the process X. In other words, for each candidate choice we simulate 
trajectories of X and consider whether, in our opinion, they provide a reasonable 
representation of our uncertainty about the evolution of the climate state. 

5 Bayesian calibration of SM91 and prediction using the 
particle filter 

We now come back to system ([5]) and operationalize the considerations of section |4] 
to determine model parameters and provide predictions. Let us observe that certain 
parameters of SM91 may be estimated by means of large numerical models of the 
atmosphere, ocean and ice sheets. We conserve the values proposed by [55 and 
ignore the time being the uncertainty on these parameters. The ai, bi and ci are 
far more difficult to estimate with general circulation model experiments and the 
sole prior constrain we have on these parameters is that they are positive. However, 
most instantiations of these parameters may be rejected by reference to the spectrum 
estimates of palaeoclimate data. Indeed, one may observe that ^ has three fixed- 
points. Keeping in mind that the longest dominant period of the forcing is 41,000 
years, the model has no chance of generating 100,000-year long cycles as those seen 
in the data, as long as the model response is quasi- linear around a fixed-point. We 
therefore need a mechanism to eject phase-space trajectories out of the fixed points 
in order to get them to explore the entire phase space. One solution is simply to 
require all three fixed points to be unstable for zero forcing |100j . 
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Figure 6: Sample paths for CO2 , 
based on choices for Tr{ijj) and S. The 
state vector was started at —450 kyr, 
using the same initial value. Panel A: 
deterministic propagation, and origi- 
nal choice of parameters ip. Panel B: 
deterministic propagation, but for 
a random sample of five different 
choices for tp. Panel C: stochastic 
propagation (five samples) , and orig- 
inal choice of ip. Panel D: Same sam- 
ple of ip as Panel B, but with stochas- 
tic propagation. 



The condition is met when the Oj, bi and Cj are set to the values published in 
|81j plus or minus roughly 20 %. We therefore use this information to build prior 
distributions for these parameters (cf. Table 1) (a more formal construction of the 
priors is planned as a next step of the present project). 

Figure |6j^ replicates the astronomically forced-evolution of CO2 starting from 
—450 kyr using the tectonic mean as initial conditions. Panel B shows five model 
realisations, starting from the same initial conditions but with different parameters 
sampled within the prior distribution. All paths show a typical oscillation of the 
order of 100-kyr — characterised by a hysteresis in the Ice vs CO2 space (Figure [7| — 
but they quickly split into groups that are out of phase with each other. The rea- 
son for this behaviour is to be found in the ability of the system to synchronise 
on the astronomical forcing, like heart beats with the pulses of a pace-maker. The 
synchronisation process is here not trivial because the astronomical forcing is aperi- 
odic (compare the autonomous and forced phase-space trajectories on Figure [7|. In 
particular, small changes in model parameters may cause synchronisation slips near 
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Figure 7: Phase-space diagrams of trajectories simulated with the SM91 model from 
—450 kyr to kyr, using the standard choice of parameters ip. The start and end of the 
trajectories are denoted by the black and white circles, respectively. The model exhibits 
a limit cycle in absence of external forcing, with a trajectory that resembles that obtained 
with data (Figure [S]) . The astronomical forcing adds a number of degrees of freedom that 
complicates the appearance of the phase diagram. The trajectory now crosses a number of 
local instability points around which a small stochastic forcing may have large effects. 



locally unstable points. This explains why originally close trajectories may become 
out of phase. 

An extra layer of complication is added by the stochastic perturbation of the 
deterministic model. Consistently with the formalism developed in section |4] we 
consider here additive perturbations of variance S. 

Figure [6p shows five sample model trajectories with stochastic forcing, but with 
the original choice of model parameters. Again, trajectories split along different 
out-of-phase groups, but compared to Figure |6p they are also less smooth. The 
combination of both parameter uncertainty and stochastic propagation is shown on 
panel D. 

The above clearly shows that sequential state and parameter data-estimation 
techniques will be needed to provide predictions in the time domain, like, for exam- 
ple, the timing of the next glacial inception. The complications introduced by the 
highly non-linear character of the model dynamics and our uncertainty on param- 
eters are such that there currently is no practical algorithm to solve this problem 
(this is a very active area of research in Statistics). Here, we used the approxima- 
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Table 1: Parameters of the physical model, choices for the stochastic components, and the 
simulation method. 

1. Fixed (non-stochastic) model parameters: 





11/253 


°C/ppmv 


kR 


0.08 


10i8°c/(Wm-2) 


ke 


0.5 


°c/°c 


Ki 


10-^ 


1/yr 


Ke 


2.5 X 10-^ 


1/yr 



2. Stochastic model parameters {iji), central values: 



ai 


8.7 X 10"^ 


10iSkg/°C yr 


hi 


1.3 X 10"^ 


1/yr 


b2 


1.1 X 10-6 


l/(ppm yr) 




3.6 X 10-8 


l/(ppm2 yr) 


he 


5.6 X 10-3 


ppm/(°C yr) 


Cl 


1.2 X 10-5 


°C/(10i8kg yr) 



For 7r(V'), the components of ij) are treated as independent, and each component 
assigned a log-normal distribution with mean equal to the central value and a 
coefficient of variation of 20% (i.e. the standard deviation is two tenth of the 
mean) . 



3. Stochastic propagation variance (S) for the state vector X = 
(/lO^8kg,/xppm,0°C): 



/0.52 \ 



S = 100-2 











V 



0.5^ 



4. Simulation method. The stochastic differential equation was simulated using 
the stochastic Euler method, using a time-step of 100 years. 
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tion proposed by Liu and West |102j based on a particle filter. It works as follows: 
At each time step, the algorithm considers n (here : n = 50000) samples of model 
parameters and states called 'particles'. Particles are propagated forward in time 
according to a three-step process: (1) forward integration of the stochastic equa- 
tions; (2) particle weight estimation as the product of prior weight and a likelihood 
taking into account the observation and its uncertainty; and (3) particle resampling 
in order keep a set of particles that all have approximately the same weight (impor- 
tance resampling). We emphasise that this algorithm is a filter and not a smoother. 
Consequently, posterior estimates at any time t are only informed by data prior to 
time t. 

Figure [8] shows one realisation of the filter. The model is initialised at —450 kyr. 
The proxy (Z) is the CO2 measured in ice-core air bubbles at Vostok, Antarctica [53] , 
but for updating we only use data between —400 and —126 kyr. CO2 presents the 
advantage of being included in the climate state vector. We simply need to specify 
an observation error, here assumed to be Gaussian with 20 ppmv standard deviation. 
Furthermore, the GT4 time scale covering the interval —400 kyr to present is only 
weakly constrained by orbital assumptions, with only two tight-points (—110 kyr 
and —390 kyr). 

This first experiment allows us to validate both the model and the particle fil- 
ter. It constitutes an interesting test to the model because —126 kyr is the end 
of the penultimate deglaciation, which puts us in conditions analogous to the early 
Holocene. First, the model correctly captures ice volume variations, here indicated 
by the SPECMAP stack record [24^ (see the Figure caption for details). The way 
the model tracks ice volume data is surprisingly good. It is not perfect, but the 
SPECMAP stack is not a perfect proxy for ice volume either because it is partly 
affected by ocean temperature. For example, the ice volume maximum at —230 kyr 
(this is event 7.4 in the SPECMAP nomenclature) is overestimated by the model, 
but there is also evidence that the event is underestimated by SPECMAP }26| I103j . 
Second, model predictions for climate evolution after —126 kyr are overall consis- 
tent with data. These predictions are shown by the thin gray lines. They are 
obtained as model realisations sampled from the latest state and parameter esti- 
mates at —126 kyr. In particular the model succeeds in predicting the draw-down in 
CO2 and the large ice mass increase around —70 kyr. Model predictions then split 
but 60 % are consistent with ice volume above 20 10^^ kg of ice at the last glacial 
maximum, and 36 % are consistent with the more stringent prediction of ice volume 
above 30 10^*^ kg of ice at that time. By contrast, the late decrease in CO2 around 
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Figure 8: Stochastic propagation of the SM91 model (initiahsed at —450, 000 years with 
priors as defined in the text), constrained by ice-core CO2 data until —126 kyr [53]. The 
training dataset is represented by triangles. Circles are validation data not used for model 
training. This includes (top) SPECMAP foraminifera data scaled as in Figure |4] and CO2 
data after —126 kyr. The filtering distributions of the state through the training period 
are represented by shades, dark and light gray representing the [25*^; 75**^] and [5*'^;95"^] 
quantiles of the particle weighted distributions, respectively. Thin lines after —126 kyr years 
represent 25 simulated predictive trajectories, started at —126 kyr using initial conditions 
and parameters sampled from the posterior estimates at that time. 
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Time 

Figure 9: Same as figure [sjbut with training until —8 kyr. Data after —8 kyr are conserva- 
tively neglected here as a test of Ruddiman's hypothesis, according to which humans have 
been perturbing climate since that time. Predicted trajectories remain, though, consistent 
with a long interglacial with no glacial inception before about 40 to 50 kyr from now, and 
a maximum in ice volume at + 60 kyr. 

—90 kyr as compared to the data constitutes a severe weakness. The relatively weak 
predictability of CO2 is also noteworthy. Poor predictability is the price we pay for 
correctly accounting for both model and data uncertainties. 

To the extent that Figure |8] seems to us to be a reasonable description of our 
uncertainty about the dynamics of Quaternary climates, we regard the combination 
of our physical model and our statistical choices as 'valid'. We are now in a position 
to provide a prediction for climate evolution after the early Holocene. Conservatively, 
any data after —8 kyr years is ignored since they may have been contaminated by 
human influence according to the early anthropogenic hypothesis. The forecast 
therefore starts at —8 kyr. 
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Results are shown on Figure [9j Again, the model provides a particularly good 
tracking of the ice volume record up to —8 kyr even though this record has not 
been used to train the model. The model prediction after the early Holocene is 
then a gradual melt-down of ice, further increase in CO2 and then a slow decrease 
up to +40 kyr, at which time glacial inception is expected with good confidence. 
Predictability after that point is poor. This prediction is consistent with the Berger 
and Loutre one [8] . 

At this stage of our research project, we still need to emphasise that the CO2 
dynamics during interglacial periods are not satisfactorily reproduced by our model, 
as we have seen on Figure |8j This implies that we cannot yet be confident that this 
prediction of the glacial inception is well-calibrated, but given the model success 
in predicting the timing of peak glacial conditions, our prediction for ice volume 
maximum at +60 kyr is, we think, trustworthy (ignoring the possible effects of an 
anthropogenic perturbation). As a straightforward extension of the present work it 
will also be informative to examine the period —800 kyr to —400 kyr, during which 
CO2 variations exhibited a smaller amplitude. It is expected that the calibration 
parameters needed to accommodate this earlier part of the record will be different, 
and it would be insightful to interpret these differences. 

As pointed out by one referee (P. Huybers), an issue may arise if structurally 
different models equally pass the validation test but yield incompatible predictions. 
In this case one must hope that an extra validation criteria will lead us to prefer one 
model to the other. The contrary will be an indication that model discrepancy was 
underestimated . 

6 Conclusion 

Behind this paper is the message that climate modelling is not and should not be 
a merely technological question, to be solved with larger and more complicated 
climate models, and faster computers. Of course, such models skillfully predict 
many complicated aspects of atmosphere and ocean dynamics; in that sense they 
are important and useful. Yet they represent just the far end of a spectrum of 
models, all of which have a role to play. 

We have emphasised the use of phenomenological models to capture the wealth of 
the data in the palaeoclimate record. These models are tractable, and a rich source 
of insights about climate behaviour on large spatial and temporal scales. In their 
complexity, they are also beautiful mathematical objects. But they present some 
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serious challenges if they are to be used quantitatively, since, with their reduced 
physics and their under-determined parameters, they make far greater demands on 
our inferential framework. We have used Bayesian statistics to combine informa- 
tion from the observational record and the model behaviour, accounting for both 
parametric uncertainty and model limitations. 

More generally, we have stressed again that the analysis of palaeo climates is 
a truly multi-disciplinary experience, starting with the field scientists who retrieve 
and interpret the observational data, and then encompassing a wide range of applied 
mathematicians, including climate modellers, complex systems experts, and statis- 
ticians. As we hope we have shown, palaeoclimates pose challenges in all of these 
areas, and questions as deceptively simple as 'When do we expect the next glacial 
inception?' require us to operate at the frontier of our knowledge and technique. 
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